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Abstract 

We develop a first-principles scheme to study ferroelectric phase transitions 
for perovskite compounds. We obtain an effective Hamiltonian which is fully 
specified by first-principles ultra-soft pseudopotential calculations. This ap- 
proach is applied to BaTiOa, and the resulting Hamiltonian is studied using 
Monte Carlo simulations. The calculated phase sequence, transition temper- 
atures, latent heats, and spontaneous polarizations are all in good agreement 
with experiment. The order-disorder vs. displacive character of the transitions 
and the roles played by different interactions are discussed. 
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Because of their simple crystal structure, the cubic ferroelectric perovskites present a 
special opportunity for the development of a detailed theoretical understanding of the ferro- 
electric phase transition. However, even in BaTiOs, a much-studied prototypical example of 
this class of compounds ||1|], many aspects of the phase behavior are far from simple. BaTiOs 
undergoes a succession of phase transitions, from the high-temperature high-symmetry cu- 
bic perovskite phase (Fig. ||) to slightly distorted ferroelectric structures with tetragonal, 
orthorhombic and rhombohedral symmetry. There is increasing evidence that the cubic- 
tetragonal transition, at first thought to be of the simple displacive kind, may instead be 
better described as of the order-disorder type. 

A comparison with the related cubic perovskites indicates that this and other aspects of 
the phase transformation behavior in BaTiOs are not universal, but rather must depend on 
details of the chemistry and structural energetics of the particular compound. Therefore, it 
is of the first importance to develop a microscopic theory of the materials properties which 
determine the ordering of the phases, the character and thermodynamic order of the transi- 
tions, and the transition temperatures. The value of a microscopic approach has long been 
appreciated, but its realization was hindered by the difficulty of determining microscopic 
parameters for individual compounds. The forms of phenomenological model Hamiltonians 
were limited by the available experimental data, leading to oversimplification and am- 
biguities in interpretation. For the perovskite oxides, empirical and nonempirical pair 
potential methods did not offer the high accuracy needed for the construction of re- 
alistic models. Recently, high quality first-principles calculations within the local density 
approximation (LDA) have been shown to provide accurate total-energy surfaces for per- 
ovskites |]7HlO[l. While an ab-initio molecular-dynamics simulation of the structural phase 
transition is not computationally feasible at present, the application of these first-principles 
methods can clearly form a foundation for the realistic study of the finite-temperature phase 
transitions. 

In this paper, we pursue a completely first-principles approach to study the ferroelec- 
tric phase transitions in BaTiOa. In particular, we (i) construct an effective Hamiltonian to 
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describe the important degrees of freedom of the system [0,|T^], (ii) determine all the param- 



eters of this effective Hamiltonian from high-accuracy ab-initio LDA calculations P, pr3| , P 
and (iii) carry out Monte Carlo (MC) simulations to determine the phase transformation 
behavior of the resulting system. We find the correct succession of phases, with transition 
temperatures and spontaneous polarization in reasonable agreement with experiment. Strain 
coupling is found to be crucial in producing the correct succession of low-symmetry phases. 
Finally, by analyzing the local distortions and phonon softening, we find the cubic-tetragonal 
transition in BaTiOs to be intermediate between the displacive and order-disorder limits. 

Briefly, the effective Hamiltonian is constructed as follows. Since the ferroelectric tran- 
sition involves only small structural distortions, we represent the energy surface by a Taylor 
expansion around the high-symmetry cubic perovskite structure, including fourth-order an- 
harmonic terms. Because the contribution to the partition function decays exponentially 
with increasing energy, we simplify this expansion by including only low energy distortions. 
Among all the possible phonon excitations, the long-wavelength acoustic modes (strain) and 
lowest transverse-optical phonon modes (soft modes) have the lowest energy. It is therefore 
our approximation to include only these two kinds of phonon excitations, thus reducing the 
number of degrees of freedom per unit cell from fifteen to six. This approximation could 
later be systematically improved, or entirely removed, by including higher-energy phonons. 

It is straightforward to describe the strain degrees of freedom associated with the acous- 
tic modes in terms of displacement vectors v/ associated with each unit cell /. In a similar 
manner, we introduce variables u; to describe the amplitude of a "local mode" associated 
with cell /. The properly chosen local mode should reproduce the soft-mode phonon dis- 
persion relation throughout the Brillouin zone, preserve the symmetry of the crystal, and 
minimize interactions between adjacent local modes. The local mode chosen for BaTiOs 
is shown in Fig. |1|. The terms in our Taylor expansion of the energy in the variables {u} 
and {v} are organized as follows: (i) a soft-mode self-energy i?'^'^'^({u}) containing intrasite 
interactions to quartic anharmonic order; (ii) a long-range dipole-dipole coupling E'^^^{{u.}) 
and a short-range (up to third neighbor) correction £''^^°''*({u}) to the intersite coupling, 



both at harmonic order; (in) a harmonic elastic energy E'^^'^^{{'v}); and (iv) an anharmonic 
strain-soft-mode couphng i?™*({u}, {v}) containing Gruneisen-type interactions (i.e., hnear 
in strain and quadratic in soft-mode variables). The cubic symmetry greatly reduces the 
number of expansion coefficients needed. All the expansion parameters are determined from 
highly-accurate first-principles LDA calculations applied to supercells containing up to four 
primitive cells (20 atoms). The calculation of the needed microscopic parameters within 
LDA for BaTiOa has been made possible by the use of Vanderbilt ultra-soft pseudopoten- 
tials |jl3| , which make large-scale calculations tractable at the high level of accuracy needed, 



and by the recent theory of polarization of King-Smith and Vanderbilt [T^, which provides 



a convenient method of calculating the dipolar interaction strengths [jT4[. The details of the 
Hamiltonian, the first-principles calculations, and the values of the expansion parameters 
will be reported elsewhere . 



We solve the Hamiltonian using Metropolis Monte Carlo simulations JT^jITII on an L x 
L X L cubic lattice with periodic boundary conditions. The homogeneous part of the strain 
in the system is separated out and treated as six extra degrees of freedom. Since most 
energy contributions (except £'*^p') are local, we choose the single-fiip algorithm and define 
one Monte Carlo sweep (MCS) as fiip attempts. 

The ferroelectric phase transition is very sensitive to hydrostatic pressure, or equivalently, 
to lattice constant. The LDA-calculated lattice constants are typically 1% too small, and 
even this small error can lead to large errors in the zero-pressure transition temperatures. 
The effect of this systematic error can largely be compensated by exerting a negative pressure 
that expands the lattice constant to the experimental value. For BaTiOs, we choose P = 
—4.8 GPa which gives the best overall agreement for the computed volumes for the four 
phases with their experimental values. The following simulations and analysis are for this 
pressure. 

In our simulation, we concentrate on identifying the succession of different phases, deter- 
mining the phase transition temperatures, and extracting qualitative features of the tran- 
sitions. We also focus on identifying the features of the Hamiltonian which most strongly 



affect tlie transition properties. For tliese purposes, it is most convenient to monitor directly 
tlie beliavior of tlie order parameter. In tfie case of tlie ferroelectric phase transition, this 
is just the polarization vector (or equivalently, the soft-mode amplitude vector u) averaged 
over the simulation cell. To avoid effects of possible rotation of the polarization vector and 
to identify the different phases clearly, we choose to accumulate the absolute values of the 
largest, middle, and smallest components of the averaged local-mode vector for each step, 
denoted by ui, U2, and u^, respectively {ui > U2 > Ms). The cubic (C), tetragonal (T), 
orthorhombic (O), and rhombohedral (R) phases are then characterized by zero, one, two, 
and three non-zero order-parameter components, respectively. As a reference, the average 
local mode amplitude u = J2i\ui\/N is also monitored. Here, Uj is the local mode vector at 
site i and is the total number of sites. 

Fig. 1^ shows the quantities ui, U2, M3, and u as functions of temperature in a typical 
simulation for an L = 12 lattice. For clarity, we show only the cooling down process. 
The values are averaged over 7000 MCS's after the system reaches equilibrium, so that the 
typical fluctuation of order parameter components is less than 10%. We find that Ui, U2, 
and M3 are all very close to zero at high temperature. As the system cools down past 340K, 
Ui increases and becomes significantly larger than U2 or M3. This indicates the transition 
to the tetragonal phase. The homogeneous-strain variables confirm that the shape of the 
simulation cell becomes tetragonal. Two other phase transitions occur as the temperature 
is reduced further. The T-0 transition occurs at 255K (sudden increase of U2) and the 0-R 
transition occurs at 210K (sudden increase of ^3). The shape of the simulation cell also 
shows the expected changes. The sequence of transitions exhibited by the simulation is the 
same as that observed experimentally. 

The transition temperatures are located by careful cooling and heating sequences. We 
start our simulation at a high temperature and equilibrate in the cubic phase. The temper- 
ature is then reduced in small steps. At each temperature, the system is allowed to relax for 
10,000 MCS (increased to 25,000 and then to 40,000 MCS's close to the transition). After 
each transition is complete, the system is reheated slowly to detect any possible hysteresis. 
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The calculated transition temperatures are shown in Table |. Simulations for three lattice 
sizes are performed; the error estimates in the table reflect the hysteretic difference between 
cooling and heating, which persists even after significant increase of the simulation time. 
(The C-T transition temperature for L = 10 is difficult to identify because of large fluc- 
tuations between phases.) The calculated transition temperatures are well converged with 
respect to system size, and are in good agreement with experiment. The saturated spon- 
taneous polarization Pg in different phases can be calculated from the average local mode 
variable. The results are also shown in Table |l[ We find almost no finite-size effect, and the 
agreement with experiment is very good for O and T phases. The disagreement for the R 
phase may be due in part to twinning effects in the experimental sample |[T9[| . 

One way to determine the order of the transition is to calculate the latent heat. An 
accurate determination of the latent heat would require considerable effort; here, we only 
try to provide good estimates. We approach the transition from both high-temperature and 
low-temperature sides until the point is reached where both phases appear equally stable. 
The difference of the average total energy is then the latent heat . This estimate should 



be good as long as some hysteresis is present. The calculated latent heat (Table |) is in 
rough agreement with the rather scattered experiment data. The discrepancy for the C- 
T transition can be partly attributed to the finite-size effect. We find that, taking into 
account finite-size effects, the latent heats for all three transitions are significantly non-zero, 
suggesting all transitions are first-order. For the T-0 and 0-R transitions, this is consistent 
with Landau theory, which requires a transition to be first-order when the subgroup relation 
does not hold between the symmetry groups below and above T^. For the C-T transition, 
although it has the largest latent heat, we find indications that it is the most weakly first 
order. The relatively large finite-size effect suggests a relatively long correlation length near 
the transition, and the change of the order parameter during the transition is also more 
gradual (Fig.^. 

Next, we investigate the extent to which the cubic-tetragonal transition can be charac- 
terized as order-disorder or displacive. In real space, these possibilities can be distinguished 
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by inspecting the distribution of the local-mode vector Uj in the cubic phase just above the 
transition. A displacive (microscopically nonpolar) or order-disorder (microscopically polar) 
transition should be characterized by a single-peak or double-peak structure, respectively. 
The distribution of Ux sXT = 4007^^ is shown in Fig. ^. It exhibits a rather weak tendency to 
a double-peaked structure, indicating a transition which has some degree of order- disorder 
character. We also see indications of this in the u — T relation in Fig.^; even in the cubic 
phase, the magnitude of the local-mode vector u is significantly non-zero and close to that of 
the rhombohedral phase. Although the components of the local modes change dramatically 
during the phase transition, u only changes slightly. 

In reciprocal space, a system close to a displacive transition should show large and 
strongly temperature-dependent fluctuations of certain phonons (soft modes) confined to a 
small portion of the Brillouin zone (BZ). For an extreme order-disorder transition, on the 
other hand, one expects the fluctuations to be distributed over the whole BZ. For BaTiOa, 
we calculated the average Fourier modulus of the soft TO mode < |M(q')P > at several 
temperatures just above the C-T transition. A strong increase of < |'u(q')P > as T — 
would indicate phonon softening. As expected, we do observe this behavior for modes at 
F. While these modes become "hard" rather quickly along most directions away from F, 
they remain soft at least half-way to the BZ boundaries along the {100} directions, again 
indicating some order- disorder character. Thus, from the example of BaTiOs, it seems that 
a positive on-site quadratic coefficient does not automatically imply a displacive character 
for the transition. Rather, the relevant criterion is the extent to which the unstable phonons 
extend throughout the BZ. 

Our theoretical approach allows us to investigate the roles played by different types of 
interaction in the phase transition. First, we study the effect of strain. Recall that the 
strain degrees of freedom were separated into local and homogeneous parts, representing 
finite- and infinite-wavelength acoustic modes, respectively. Both parts were included in the 
simulations. If we eliminate the local strain (while still allowing homogeneous strain), we 
find almost no change in the transition temperatures. This indicates that the effect of the 
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short- wavelength acoustic modes may not be important for the ferroelectric phase transition. 
If the homogeneous strain is frozen, however, we find a direct cubic-rhombohedral phase 
transition, instead of the correct series of three transitions. This demonstrates the important 
role of homogeneous strain. Second, we studied the significance of the long-range Coulomb 
interaction in the simulation. To do this, we changed the effective charge of the local 
mode (and thus the dipole-dipole interaction), while modifying other parameters so that the 
frequencies of the zone-center and zone-boundary phonons remain in agreement with the 
LDA values. We found only a slight change (10%) of the transition temperatures when the 
dipole-dipole interaction strength was reduced by half, but a further reduction changed the 
results dramatically (in fact the ground state becomes a complex antiferroelectric structure). 
This result shows that it is essential to include the long-range interaction, although small 
inaccuracies in the calculated values of the effective charges or dielectric constants may not 
be very critical. On the other hand, our tests do indicate a strong sensitivity of the Tc's to 
any deviation of the fitted zone-center or zone-boundary phonon frequencies away from the 
LDA results. Thus, highly accurate LDA calculations do appear to be a prerequisite for an 
accurate determination of the transition temperatures. 

Our approach opens several avenues for future study. Allowing a higher-order expansion 
of the energy surface might allow an accurate determination of the phase diagram. More 
extensive Monte-Carlo simulations on larger systems, and with careful analysis of finite-size 
scaling, could provide more precise transition temperatures, free energies, and latent heats 



| 21| . Finally, the theory would be more satisfying if the 1% underestimate of the lattice 
constant in the LDA calculation could be reduced or eliminated. 

In conclusion, we have developed a first-principles approach to the study of structural 
phase transitions and the calculation of transition temperatures in BaTiOs. We have ob- 
tained the transition sequence, transition temperatures, and spontaneous polarizations, and 
found them to be in good agreement with experiment. We find that long- wavelength acous- 
tic modes and long-range dipolar interactions both play an important role in the phase 
transition, while short-wavelength acoustic modes are not as relevant. The C-T phase tran- 
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sition is not found to be well described as a simple displacive transition; on the contrary, if 
anything it has more order-disorder character. 

We would like to thank R.D. King-Smith, U. V. Waghmare, R. Resta, Z. Cai, and 
A.M. Ferrenberg for useful discussions. This work was supported by the Office of Naval 
Research under contract numbers N00014-91-J-1184 and N00014-91-J-1247. 
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TABLES 

TABLE I. Calculated transition temperatures Tc, saturated spontaneous polarization Pg, and 

estimated latent heat /, as a function of simulation cell size. 
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FIGURES 

FIG. 1. The structure of cubic perovskite compounds BaTiOs. Atoms Ba, Ti and O are 

represented by shaded, solid, and empty circles respectively. The areas of the vectors indicate the 
magnitudes of the displacements for a local mode polarized along x. 

FIG. 2. The averaged largest, middle, and smallest components ui, U2, U3 and amplitude u of 
local modes as a function of temperature in a cooling-down simulation of a 12 x 12 x 12 lattice. 
The dotted lines are guides to the eyes. 

FIG. 3. The distribution of a Cartesian component of the local mode variable in the cubic 
phase at T = 400K. 
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